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We present a detailed analysis of the expected signal-to-noise ratios of supermassive black hole 
binaries on eccentric orbits observed by pulsar timing arrays. We derive several analytical rela¬ 
tions that extend the results of Peters and Mathews |l| to quantify the impact of eccentricity in 
the detection of single resolvable binaries in the pulsar timing array band. We present ready-to-use 
expressions to compute the increase/loss in signal-to-noise ratio of eccentric single resolvable sources 
whose dominant harmonic is located in the low/high frequency sensitivity regime of pulsar timing 
arrays. Building upon the work of Phinney @| and Enoki and Nagashima Qj, we present an ana¬ 
lytical framework that enables the construction of rapid spectra for a stochastic gravitational wave 
background generated by a cosmological population of eccentric sources. We confirm previous find¬ 
ings which indicate that, relative to a population of quasi-circular binaries, the strain of a stochastic, 
isotropic gravitational wave background generated by a cosmological population of eccentric binaries 
will be suppressed in the frequency band of pulsar timing arrays. We quantify this effect in terms 
of signal to noise ratios in a pulsar timing array. 

I. INTRODUCTION 

It is believed that supermassive black holes (SMBHs) 
with masses between 10 6 Mq — 1O 9 M 0 are ubiquitous in 
galactic nuclei M- According to the accepted frame¬ 
work of hierarchical structure formation, massive galax¬ 
ies are formed by continuous accretion of gas from cosmic 
web filaments or through galactic mergers This lat¬ 

ter mechanism naturally leads to the formation of SMBH 
binaries in the merged galaxy remnants. As the SMBHs 
sink in the potential well of the remnant galaxy due to dy¬ 
namical friction, stars within the binary orbit are quickly 
ejected. An SMBH merger can only take place if addi¬ 
tional mechanisms operate to remove energy and angular 
momentum from the binary, e.g., friction from a spheri¬ 
cal Bondi accretion flow [f|, a circumnuclear gas disk 0, 
slingshot scattering of stars on low angular momentum 
orbits intersecting the binary EMI , etc. If any of these 
mechanisms can drive the orbit to sufficiently small sepa¬ 
rations, gravitational wave (GW) emission can take over 
and drive the binary system the rest of the way to coa¬ 
lescence within a Hubble time El: H~il HU . 

Regarding the orbital properties of SMBH binaries, 
scattering interactions between individual stars and 
SMBH binaries can potentially drive the binaries to 
large orbital eccentricities, particularly when the bina¬ 
ries retain significant eccentricities at the end of the dy¬ 
namical friction phase El, mUIll, whereas SMBH bi¬ 
naries embedded in sufficiently massive prograde self- 
gravitating gas disks may acquire eccentricities as large 
as e ~ 0.6 — 0.8 by the time gravitational radiation takes 
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over the dynamical evolution of the system E3- Fur¬ 
thermore, SMBH binaries embedded in counter rotating 
disks may be driven to very large values of eccentricity 
e ~ 1 [23, |23| , even though the binary can flip and realign 
with the disk [24j . 

The gravitational radiation emitted during the inspi¬ 
ral of binaries with masses 1O 6 M 0 — 1O 9 M 0 out to red- 
shifts z <1 will be detectable by Pulsar Timing Arrays 
(PTAs) l25l - l34l |. PTAs are capable of detecting cosmic 
string networks, primordial GWs, an unresolved stochas¬ 
tic GW background generated by a large population of 
compact binary sources !~iT)l| and GWs from individual 
binary systems j4l|, [HI . 

Given the significant attention that eccentric com¬ 
pact binaries have attracted as potential sources of GWs 
and electromagnetic radiation |43H47T | , there is a need to 
study the effect of eccentricity both in terms of source 
detection and parameter estimation for individually re¬ 
solvable sources, and for the detection of a stochastic GW 
background in the context of PTAs. Our understanding 
on the effect of eccentricity on potential GW sources for 
PTAs has gradually improved from the seminal work of 
Quinlan |l2|, and recent theoretical and numerical stud¬ 
ies that have shed light on the impact of eccentricity and 
environmental effects in suppressing the low frequency 
GW background in the PTA band [48M55j . 

In this article we build upon the work of Phinney Q 
and Enoki and Nagashima [ 3 ] by constructing an ana¬ 
lytical framework that enables the construction of rapid 
spectra for a stochastic GW background generated by a 
population of eccentric sources. We then employ a pre¬ 
scription for the evolution of the BH mass function taken 
from [56|, and combine it with our results to compute the 
signal-to-noise ratios (SNRs) of a stochastic GW back¬ 
ground generated by a population of eccentric binaries, 
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with the SNR in that case being derived from a cross¬ 
correlation statistic, given that matched filtering cannot 
be applied to a stochastic signal. We also derive several 
analytical summations that expand upon the results of 
Peters and Mathews [l] (PM hereafter) to explore in de¬ 
tail the effect of eccentricity on the GW strain and the 
matched-filter SNRs of individually resolvable sources. 

Our studies show conclusively that the SNR of eccen¬ 
tric binaries is non-negligibly attenuated for eccentricity 
values (e ;> 0.7). However, binaries with low to moderate 
values of eccentricity (0 <; e <; 0.6) will have SNRs com¬ 
parable to their quasi-circular counterparts. This sug¬ 
gests, in principle, that the detection of a population of 
eccentric binaries may be possible and would provide new 
insights on the formation channels of SMBH binaries and 
their cosmological evolution. However, it is still necessary 
to show that the imprints of eccentricity can be accu¬ 
rately extracted from GW observations with PTAs. We 
defer the study of this important issue to future work. 

This article is organized as follows: in Section |TT| we 
provide a succinct description of the properties of eccen¬ 
tric binary systems and derive analytical relations that 
are of importance for eccentric SMBH binaries observed 
by PTAs. In Section [ITT] we provide analytical results 
for the energy density and the characteristic amplitude 
of the GW spectrum, and discuss at length the effect of 
eccentricity on these two observables. In Section [IV] we 
apply this calculated strain to compute SNRs for both 
single resolvable sources and a stochastic population of 
eccentric binaries with e € [0, 0.9]. We summarize our 
findings and describe future directions for research in Sec¬ 
tion |Vj Throughout this article we use geometric units 
with G = c = 1. 


II. POWER FROM INDIVIDUAL ECCENTRIC 
BINARIES 


Consider a binary system with component masses 
(mi, m 2 ), such that mi > m 2 , M = mi+m 2 , and whose 
orbital rest-frame frequency is given by / orb = w/27t. If 
the system evolves from an initial state with nonnegligi- 
ble eccentricity e and semi-major axis a, then the binary 
radiates GWs in the whole spectrum of harmonics. Fur¬ 
thermore, as shown by PM [l|, the relative power P(n) 
radiated in the n’th harmonic is given by: 


p (n) 


32 mfm| (mi + m 2 ) 
5 a 5 


9 {n,e ), 


(1) 


g(n, e) = 


32 


Jn— 2 (^.e) 2eJ„_i(ne) 


( 2 ) 


T Jnij ^e) T 2eT^_(_i(ne) Ty l _(_2(ne) 


+ (l - e 2 ) | J„_ 2 (ne) - 2J„(ne) + J„ +2 (ne)| 


3n 2 


J 2 (ne) 


Using Bessel’s equation and recurrence relations, one can 
re-write Eq. © as follows: 



Note that Eq. © corrects a typo in Eq. (Al) of PM [l|. 
As shown in PM: 


F ( e ) = Y e ) 

n—1 


1 + —e 2 4- — e 4 

' 24 ' 96 

(1 - e 2 ) 7/2 


( 4 ) 


Hence, averaging over one period of the elliptical motion, 
the average rate at which the binary system radiates en¬ 
ergy is given by: 


OO 

<P> = £P(n), 

n=1 

32 m?TOo M / 73 2 37 4 \ 

5 a 5 (! — e 2 ) 7/2 ( + 24® + 96 6 ) ' 


( 5 ) 


Another interesting quantity that involves the object 
g(n , e) is the GW strain root-mean-square (rms) ampli¬ 
tude. As discussed in [57], the rms amplitude and the 
energy radiated in the n’th harmonic are related through: 


1 + z \] E n 

Trdl, Tl /orb 


( 6 ) 


where z is the redshift. Since the luminosity, E, emitted 
by the system averaged over one complete orbit is given 
by 


09 00 

E = -^-M 10/3 (27r/ orb ) 10/3 9 (n, e ), (7) 


where 
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then Eq. © can be re-written as follows: 

Io o a^ 5/3 

hn = 2 W— - — (27r/ orb ) 2/3 \/g( n i e)(l + 2), (8) 

V 5 naL 

where M. = M rj 3 / 5 is the chirp mass, and 7? = mirri 2 /M 2 
represents the symmetric mass ratio. It is possible to ob¬ 
tain a similar expression to the average power by consid¬ 
ering the quantity 




32 M 10/3 
T d 2 


(2^/orb) 473 ^ 
n= 1 


g(n, e) 
(n/2) 2 ’ 


(9) 


where d = di/(l + z). This quantity has heretofore been 
evaluated numerically using a given number of harmonics 
to ensure a specified accuracy. However, one can derive 
an exact closed form for the sum appearing on the right- 
hand side of this expression, as shown in Appendix [A] 


H(e) = ^ 

n= 1 


g(n , e) 


4 - y/1 - e 2 
12V1 - e 2 


Thus, Eq. © takes the simple form: 


( 10 ) 


they will no longer be instantaneously monochromatic. 
Instead, we can regard the emission at each harmonic to 
represent a separate population of sources. Based on this 
observation, and following Enoki and Nagashima [3}, we 
have that: 


with 


£gw = ^ Gw > n ’ 

n— 1 

nOO pOO ^ 

£GW,n= / / N(z)——f n — , 

Jo Jo 1 + 2 d/„ )r 


d-EQw ,n | df n 

d z- 


(14) 


where f n represents the frequency of the n’th harmonic 
observed today, and /„ ir = (1 + z)f n is the frequency of 
the harmonic in the rest frame. The amount of energy 
radiated in GWs into the n’th harmonic as the frequency 
of the n’th harmonic changes from /„ >r to f n ,r + d f n>r is 
given by: 


dEcw.n 

d fn,r 


d/n,r • 


(15) 
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In the following Section, we will use a similar approach to 
derive new analytical relations to explore the signatures 
that a population of eccentric binaries may imprint on a 
stochastic background of gravitational radiation and on 
single resolvable sources. 


III. STOCHASTIC BACKGROUND OF A 
POPULATION OF ECCENTRIC BINARIES 


Following Ref. @|, one can define the total GW energy 
density per logarithmic frequency interval observed today 
from a population of (instantaneously monochromatic) 
sources as: 


£gw = 




(/) 


N(z) 


o 

1 , dEcw 

1 + z h dfr 


d/ 
/ 


d 2y. (12) 


The rate of mergers per unit comoving volume which 
occur between redshift z + dz is given by N(z)dz. Fur¬ 
thermore, p c represents the rest-mass energy that would 
be required to close the Universe [§) 


Pc = 


3Hl 

8n 


(13) 


In practice, we can replace N(z)dz by a differential rate 
and integrate over source parameters, but for the mo¬ 
ment we shall assume that the population is composed 
of identical sources. If the sources have eccentricity then 


This outgoing energy is measured in the source’s rest 
frame, and is integrated over the entire radiating lifetime 
of the source and over all solid angles jl}. Using the 
relations 

47T 2 / 2 rb a 3 = M and to = 2n f olh , (16) 

with Eq. ©, we find that: 

dEgw,n _ d^GW.n dt r 

d f n ,r dt r d f n>r ’ 

= f(A(A 0/ U ( ". e) , (18) 

ihi = jks/ 3 „h /3 (19) 

dtr dt r 5 27T 

where F(e) was defined in Eq. ©. Combining these, we 
find 


d-Eow.n = (2tt) 2/3 A4 5 / 3 g(n,e) -i 
dfn,r ~ 3 F(e) n /orb 

= ^ .(. ,e) ■ 

3F(e) (1 + z ) 1/3 (n/2) 2/3 

The energy density in the background per logarithmic 
frequency interval is then given by 

Pc^Gw(f) = jf 2 h 2 c (f) = (21) 

-Ad 5/3 (tt /) 273 [°° y, 1 g (n,e) N(z) 

3 Jo h F ^H2) 2/3 (l + z) 1/3 Z ' 

Note that in the quasi-circular limit (n = 2, e —> 0), 
Eq. (1?T1) recovers the results presented in Ref. [2]. 
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A. Estimating the number of merger events in unit 
comoving volume N(z) 


One important ingredient in the calculation of the 
stochastic spectrum, energy density and, ultimately, the 
SNR with which a population of GW sources can be de¬ 
tected is the number of mergers that occur between red- 
shift 2 and 2 + dz, i.e., N(z). For the systems under 
consideration, i.e., binaries with total masses between 
10 6 ~ 9 Mq, our knowledge of the numbers and mass dis¬ 
tributions of SMBHs has changed considerably with the 
advent of large scale surveys [HU, [59j and recent theoreti¬ 
cal studies [Hfl, inm. However, deriving a robust model 
for the computation of N(z) is a complex problem due 
to the large uncertainties inherent in several aspects of 
the calculation, e.g., the poorly constrained rate of BH 
migration toward the center of merging galaxies caused 
by interactions with dark matter, gas, and stars; the pos¬ 
sibility of multiple BH interactions in the event that the 
BH migration is inefficient, etc. [HI [5li ]62j-[73j]. With 
these caveats in mind, we use the estimate for N(z) de¬ 
scribed in Ref. [56j ] , which we will review here for com¬ 
pleteness. 

We need to estimate the comoving number density of 
BHs with masses between M. and M, + dM«. BH masses 
are strongly correlated with the bulge masses of their 
hosts and so this is equivalent to considering the distri¬ 
bution of galaxy bulge masses. The starting point for 
such an estimate is an empirical model known as the 
Schechter function [74], given by 

< j>(M)dM = ip M a exp (— M) d M , (22) 


where ip and a represent the normalization of the lumi¬ 
nosity function and the faint-end slope parameter, re¬ 
spectively. The Schechter function is a power law that is 
truncated at large masses. For the most massive galaxies 
of interest, we need to amend this function to account for 
the observed excess of mass in the brightest cluster galax¬ 
ies and other very massive elliptical galaxies. Following 
Ref. f6lj], we do this by adding a Gaussian component to 
Eq. (E2D: 


= (tp + <p 
+ <f> exp 


massive 


) d M = ip M a exp (-M) d M 
2 

' d M. 


1 / 2.5 log M\ 


(23) 


where <f> and ip are normalization factors to describe the 
brightest cluster galaxies and less massive galaxies, re¬ 
spectively. We try to encapsulate in a conservative way 
the current knowledge we have from galaxies that host 
BHs with masses ~ 10 9 Mq such as M87. Hence, follow¬ 
ing Ref. [56| we set cj> = ip and a = 0.58, which ensures at 
least one M87-mass source in our sample. The comoving 
density of BHs can be constructed from the Schechter 
function by replacing 


M 


M. 

~M~ 


with 


1.2 x 10 8 


Mr,. 


(24) 


Here M, denotes the BH mass and M is a Schechter pa¬ 
rameter that represents the characteristic mass at the 
turnover of the mass function. This particular prescrip¬ 
tion is consistent with observational data [75(. We set 
the normalization of the luminosity function to have the 
constant value 



FIG. 1 . Redshift evolution of the black hole mass function 
given by Eq. Oil. 


ip = 3 x 10 -3 Mpc -3 . (25) 


This choice is in good agreement with results presented 
in Ref. [76] at low redshifts and using the cosmological 
parameters presented in Ref. J59|. Observational data 
suggests that ip might have a mild dependence on red- 
shift. However, following Ref. [56|, we ignore the red- 
shift dependence of ip because it is a small effect that has 
a negligible influence on the total GW signal. Finally, 
ensuring that the redshift dependence of the faint-end 
slope parameter a satisfies mass conservation, one finds 
that [56| : 


a 


-2 + 


0.52 
1 + 2 ' 


(26) 


We can reconstruct the BH mass function introduced in 
Ref. [Efi] by plugging Eqs. (l24l) (l26l) into Eq. (l23l) see 
FigureQ] This approach reproduces the results presented 
in Ref. [60] at a 2cr level. Following Ref. [56|, we express 
the number density of mergers N(z) by assuming that 
it is proportional to the product of the number density 
of the constituent black holes, as shown in Eq. (8) of 
Ref. ■ Using this approach, we evaluate the integral 


N 0 = 



N(z) 

(1 + 2) 1/3 


d2, 


(27) 


where 2 m i n = 0 and 2 max = 1. Assuming that all systems 
in the Universe have the same eccentricity we find that 


f 2.63 x 10 -3 Mpc -3 , 7 <; logM <; 7.9 , 
| 1.16 x 10 -3 Mpc -3 , log M £ 7.9 . 


1 + 2 


(28) 
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We have used the mass ranges quoted above motivated by 
the fact that mass function has a break around log M ~ 
7.9 for all redshifts of interest. Note that even though 
this is a rough approximation, we have verified that this 
choice does not have a strong influence on the results. 


B. Ready to use expressions for the gravitational 
wave energy density and the characteristic 
amplitude of the gravitational wave spectrum 


Having derived all of the ingredients to compute the 
GW energy density and the characteristic amplitude of 
the GW spectrum, and using the most recent results for 
the cosmological parameters released by the Plank Col¬ 
laboration in Ref. m to compute the critical density of 
the Universe defined by Eq. m, we can derive ready- 
to-use expressions for the energy density and the char¬ 
acteristic amplitude of the GW background. We carry 
out this calculation in two steps. We first provide a ped¬ 
agogic example in which the eccentricity of the SMBH 
binary population is assumed to be constant. Thereafter, 
we address the likely physical scenario in which the ec¬ 
centricity evolves as a function of frequency due to GW 
emission. 


1. Compact binary population with fixed eccentricity 


Assuming that the eccentricity of the binary popula¬ 
tion is fixed, we can derive an analytical expression that 
reproduces the sum in Eq. m to better than 0.01% in 
the eccentricity range eo € [0, 0.95] (see Appendix ICl). 
namely: 


00 / _ \ i_i 1467 „2 _ 115 „4 , 227 „6 

A! \ 9\ n i e 0) __ 1 ' 1024^0 12288 e 0 32768^0 

(0j k{n/2f* (1-egf/ 2 


(29) 


For later convenience, let us define the function: 

A( P \ (1 i 1467 2 _ _+15_ 4 , _227_ 6\ 

TD/„ ^ _ A \ e 0) „2\ V 1 ' 1024^0 12288 e 0 3 32768 e 0 ) 

B(e »> = FM - (' - -1 + Ieg + geJ- 

(30) 

Using Eq. m, the energy density and the characteristic 
amplitude of the GW background take the form: 


ficw(f) = 3.6 x 10 


-10 


M 


10 8 Mq 


5/3 


Nn 


10- 3 Mpc -3 


33 ) B{e o), 


h c (f) = 5.0 x 10 


-16 


M 


10 8 Mq 


5/6 


No 


10-3 Mpc -3 


1/2 


lyr“ 


lyr- 1 


VB(e 0 ) ■ 


2/3 

(31) 

-2/3 

(32) 


In Figure [2] we plot the attenuation function B(e o) (see 
Eq. (1301) '). We notice that both the energy density and 
the characteristic amplitude of the GW background are 
maximized for a population of quasi-circular binaries, 
and steadily decrease for increasing values of eccentricity. 
These results give the energy density and typical strain 
of a GW background generated by binaries with fixed 
eccentricity and chirp mass. 



FIG. 2. Attenuation factor, B(e o), as defined in Eq. (1301) . 
which describes the decrease in the emitted energy density 
for a population of compact sources with fixed eccentricity. 
In light of Eqs. m and dm the present-day energy den¬ 
sity (Igw (/) is maximized for a population of quasi-circular 
compact binaries, whereas its value is decreased by a factor 
~ 10 for a population of highly eccentric systems (eo ~ 0.9). 
Similarly, the characteristic amplitude of the GW spectrum 
steadily decreases as the eccentricity of the compact binary 
population increases. 


2. Compact binary population with evolving eccentricity 

To describe a compact binary population whose ec¬ 
centricity is evolving, we notice that for a given initial 
eccentricity eo at a fiducial initial orbital frequency /o, 
the eccentricity depends only on the orbital frequency: 
e = e(/ or b, eo). Each harmonic n contributes to the sig¬ 
nal at an observed frequency / = n f OI b/(l + z). Hence, 
including the frequency evolution of the eccentricity en¬ 
tails replacing the argument of the g(n, e), F(e) functions 
in Eq. (I2T1) by 


^(/orbi eo) — e 



(33) 


We shall use the dictionary e —> e(/ or b) given by Eq. 
(3.12) of Ref. Q, which is robust for e G [0, 0.9], namely: 


^(/orbi 6 q) i 


16.83 — 3.814/3°- 3858 
16.04 + 8.1 /3 1 - 637 ’ 


(34) 
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FIG. 3. The panel shows the frequency evolution of the func¬ 
tion S(f, fo,eo,z), given by Eq. (l36l) for several values of ini¬ 
tial eccentricity eo. The ir-axis shows the ratio f r /fo, where 
f r = (1 +z)f. 


where /3 = X 2/3 / a o and x = forb/fo, with / orb = 
(1 + z) f/n, and 


cr 0 = 


12/19 

U)_ 

1 — e o 


1 + 


121 

304 ( 


870/2299 


(35) 


We substitute Eq. (1331) into Eq. m to obtain the func¬ 
tion: 


S(f,fo,e 0 = 0,z) = l, 
lim/->oo S(f, f 0 , e 0 , z) -*■ 1. 

In light of this analysis, we have constructed a function 
that has these generic properties. We found it convenient 
to split the function in two pieces given its distinct prop¬ 
erties before and after it reaches 5=1. The points at 
which S(f, f 0 ,e 0 ,z) = l are given by: 


a: fixed 


3620 e 0 370 2 i 132 

841 (1 - elf V ~ 243 6 ° + ' 


(39) 


In the domain x ;> a; fixed , we propose the following ansatz 
for S(f,f 0 ,e Q ,z): 


S high (f,fo,e 0 ,z) = l + a(e 0 )[x-b(e 0 )]e c(e ° )x , (40) 

where x = f r /fo and the eccentricity dependent coeffi¬ 
cients a(e o), 6(eo), c(eo) are given by: 

b(e 0 ) = * fixed , (41) 

C ( e °) = 2 -max _ ^fixed 1 (^) 

cm ax i 

a ( e °) = j.max _ fixed eX P ( C ( e o) *““) • (43) 

It is worth pointing out that for low values of eccentric¬ 
ity (eo <; 0.2), Eq. HU1) reproduces the main features of 
S(f, fo, eo, z) throughout the domain x > 1. When we 
consider eo > 0.2, we need to replace the low frequency 
evolution using the following relation 


S(f,fo,eo,z) 


E 


1 g(ra,e(/orb; e 0 )) 
E(e(/ OT b; e 0 )) (n/2) 2/3 


(36) 

In Figure [3] we show the frequency evolution of the func¬ 
tion S(f, fo, eo, z ) for several values of initial eccentricity 
eo- It is worth pointing out that these results are in excel¬ 
lent agreement with Ref. [3j, even though we have used a 
different approach to parameterize the orbital frequency 
evolution. We have found several interesting properties 
of the generating function S(f, fo,eo,z): 


• The location of the maxima follows a simple rela¬ 
tion given by: 


1293 

TsT 


,12/19 

'_0 _ 

1 e 0 


, 121 2 
1 + 304 e ° 


870/2299' 


3/2 


(37) 


where x = f r /fo- 

• The maxima of the S(f,fo,eo,z) function is the 
same for all values of eo and is given by 

5(/,/ 0 ,e 0 ,z) max = —. (38) 

• Two additional properties that S(f,fo,eo,z) must 
satisfy are: 


Siow(/, fo, eo, Z) = d(e 0 ) x^~ s ^ 7 e~ a ^ x , (44) 

where the coefficients (d(e o), s(eo), gieo )) are determined 
by enforcing that Slow has the correct value at x = 1 and 
x = z fixed , and that S 1 , ow (a; fixed ) = S^x^). The 
transition from S\ ow to Shigh is at the point x fixed . 

We have found that Shi g h(/, fo, e o, z), given by 
Eq. (H01) , can accurately describe the full numerical solu¬ 
tion of Eq. (l36l) for eo £ [0, 0.9] in the domain x £ x fixed . 
This is possible because the numerical solution has self¬ 
similarity properties that are captured by Eqs. Em m- 
We have attempted to provide a similar parameteriza¬ 
tion for the spectra in the domain 1 <, x <; a: fixed and 
have found that self-similarity is present for populations 
with eo 0.7. Populations with larger eccentricities 
have two properties that deviate from self-similarity in 
the domain 1 <; x <; £ fixed : (a) the slope of the spectra 
evolves as a function of eccentricity; (b) the spectra de¬ 
velops a bulging at lower frequencies that becomes more 
pronounced for increasing values of eccentricity. These 
two properties are clearly shown in the bottom panel of 
Figure |3] The parameterization we propose in Eq. HI) 
captures the evolution of the spectra as a function of ec¬ 
centricity with the parameter s(eo). Using both Eqs. HU1) 
and (1441) , we can analytically reproduce the numerical so¬ 
lution of Eq. (1331) for systems with eccentricity eo <, 0.7 
with an accuracy better than 10% in the domain x> 1 
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note that the largest deviation between the numerical and 
analytical solutions occurs for populations with eo = 0.7. 
The discrepancy arises because, even if we have captured 
the evolution of the slope of the spectra as a function 
of eccentricity, the numerical solution presents an ad¬ 
ditional bulging at low frequencies that is not equally 
present in all spectra. Indeed, the bottom panel of Fig¬ 
ure 2] shows that this feature becomes increasingly pro¬ 
nounced for highly eccentric populations in the low fre¬ 
quency domain. However, we notice that S\ ow (f, fo, eo, z) 
still provides an approximate description of the spectra 
in the low frequency domain that smoothly asymptotes 
to the numerical solution when x —> x fixed . This is an 
important property, since this is the region where the 
signal is most likely to be detected. Therefore, given the 
ever-increasing attenuation of the spectra for very large 
eccentricities, it seems that the analytical framework we 
have constructed covers the entire domain of detectable 
stochastic signals. Finally, we note that, by construction, 
our analytical approach satisfies S(f,fo,eo = 0, z) = 1 
and Imif-yoo S(f, fo,eo,z) —> 1. In future studies that 
aim at modeling SMBH binaries that evolve in stellar en¬ 
vironments or embedded in counter rotating disks that 
may drive the eccentricity to large values eo ~ 1, it will 
be necessary to modify the framework described above 
by including a non-self-similar evolution for the low fre¬ 
quency evolution part of the spectrum, in particular for 
eccentricities eo ^ 0.7. 

The analytical approximation to the spectra of eccen¬ 
tric populations we have constructed above provides a 
robust description of the imprint of eccentricity over a 
wide range of parameter space. Given its simplicity, it 
provides an ideal tool to be implemented in detection 
pipelines. We utilize this result in the following Section 
to compute the expected signal-to-noise ratio with which 
a cosmological population of SMBH binaries with non- 
negligible eccentricity can be detected with PTAs. 


IV. SIGNAL-TO-NOISE RATIOS FOR PULSAR 
TIMING ARRAYS 


In this Section we discuss in detail the prospects of de¬ 
tecting a cosmological population of inspiralling SMBH 
binaries with PTAs. Current studies suggest that the 
expected signal from these events may comprise a su¬ 
perposition of two distinct contributions: (i) a stochas¬ 
tic background generated by the incoherent superposi¬ 
tion of gravitational radiation emitted from the whole 
SMBH population [4^, [z!, Hcj; and (ii) GW signals that 
stand above the background and can be individually re¬ 
solved [H, |8I]. The motivation to consider these two 
complementary cases stems from the fact that an inho- 
mogenous combination of multiple sources emitting in 
the same frequency bin can adopt several configurations 
in the timing residuals, such as a nearly isotropic distri¬ 
bution over the sky or a few bright spots in the sky if they 
superpose coherently [45 ■ There has been a vigorous re- 
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FIG. 4. The top panel shows a direct comparison between 
the numerical solution of the sum S(f,fo,eo,z) given by 
Eq. HU and the analytical solution we have constructed using 
Eqs. (HD and (144)1. This analytical solution reproduces the 
full numerical solution for systems with eccentricity eo < 0.7 
with an accuracy better than 10% in the domain x > 1. The 
largest deviation occurs for populations with eo = 0.7, which 
start to deviate from self-similar solutions in the low frequency 
regime (1 < x < * fixed ). Bottom panel: populations with 
higher eccentricity have a non self-similar evolution in the 
domain 1 < x < * fixed , but we can still provide an approx¬ 
imate description in this regime using Eq. (ITU . Please note 
that Eq. provides a reliable description of the spectra 
for any value of eccentricity eo € [0, 0.9] for x > x flxed , and 
Eq. (144II smoothly asymptotes to the numerical solution when 
x ->■ x Rxed . 


search program to develop data analysis techniques in 
the limiting cases of an isotropic stochastic background 
which, as described in the previous Section, may be de¬ 
scribed by a power law spectrum [35 - 38 L 82l. |83| . for sin¬ 
gle monochromatic GW sources p33l [8414871 and, more re¬ 
cently, for anisotropic GW backgrounds [88l490| , although 
we will not discuss these further here. 


























A. Sensitivity of PTAs to single resolvable sources 
and a stochastic gravitational wave background 


B. Signal-to-noise ratio calculations for single 
resolvable sources 


The sensitivity curves of PTAs to continuous waves 
and a stochastic GW background have been discussed 
at length in Ref. 00. If we define er rms as the rms 
timing noise, and 1/At as the cadence of the measure¬ 
ments, then combining Eqs. (40) and (42) of Ref. 9l|, 
the dimensionless effective noise amplitude for the tim¬ 
ing residuals induced for a stochastic GW background is 
given by: 

h%(f) = fSn(f) = 24tt 2 A t a r 2 ms f . (45) 


Assuming a total observation time T 0 b s , the analysis pre¬ 
sented in [0J shows that the power law integrated sen¬ 
sitivity curve for a PTA’s response to a stochastic GW 
background has a sharp cut-off in sensitivity at a fre¬ 
quency / = T7 . On the other hand, for individually 
resolvable sources, the maximum sensitivity is attained 
around frequencies T^, and there is a slow diminish¬ 
ing in sensitivity below this value. Assuming a quadratic 
timing model, Ref. [3l| shows that the dimensionless ef¬ 
fective noise amplitude for continuous waves can be mod¬ 
eled as a two-part power law in /. This two-part power 
law, as given in Ref. [30 ] . is formally a continuous sum 
of the two components, but it will prove convenient for 
us to approximate it as a piecewise combination of the 
components, namely: 


ftc,hi g h(/) =Bf* , for 

tol 

h c ,iow(f) = c f~2 , with 


B = 
C = 


f 36 
\N P (N P 
8 B 


7~>3 
1 obs 


i) 



^rms ? 


(46) 

(47) 

(48) 

(49) 


The quantity h c (f) is the characteristic strain of noise 
fluctuations in the detector, which is required to com¬ 
pute the SNR using Eq. m below. We note that 
the transition frequency value /trans = 2 T~^ s at which 
h c , high(/) = h c ,\ow{f) is simply an approximate value 
for which the two-part power law representation of the 
total sensitivity reproduces a fully numerical Bayesian 
analysis a 

We emphasize that the effective sensitivities above dif¬ 
fer depending on the detection statistic being assumed, 
and this has occasionally resulted in some confusion 
when calculating sensitivity curves, particularly their 
spectral slopes, throughout the literature. We have as¬ 
sumed in Eq. (l45l) that the stochastic background is 
searched for using a cross-correlation statistic, whereas in 
Eqs. (H7>1) HMD, we assume that continuous-wave sources 
are searched for using matched filtering. 

Having described the prescription we will use for the 
sensitivity of PTAs to detect continuous wave sources and 
a stochastic GW background, we will now compute the 
expected SNR of single resolvable sources. 


Several recent studies have explored the ability of 
PTAs to resolve GW sources individually. For instance, 
assuming the existence of a population of quasi-circular 
monochromatic sources, an array of pulsars which are 
equally-sampled every two weeks for ten years, and mak¬ 
ing several other simplifications regarding the nature of 
the data sets, Ref. [92[ concluded that N s sufficiently loud 
sources with SNRs ;> 10 can be resolved and localized in 
the sky with a network of 3 N s pulsars. Building upon 
this study, Ref. @ demonstrated that it was possible 
to: (i) recover the SNR of injected signals to within a 
few percent; (ii) infer the sky localization to within a few 
degrees; and (iii) resolve the frequency at which the sig¬ 
nals were injected to better than 0.1 nHz. To put this 
latter result in context, a PTA that collects data for a 
time span of T 0 bs cannot in principle distinguish two GW 
frequencies separated by less than A/ ~ l/Tobs ~ 3nHz 
for T 0 b s = 10yr. If the algorithm introduced in Ref. (93) 
is capable of determining sub-Fourier bin precision to the 
level of 0.1 nHz, this means that they are capable of re¬ 
solving up to 30 sources per frequency bin. 

A more conservative approach to estimate the number 
of GW sources that can be individually resolved with a 
PTA was presented in Ref. j94j . Basic counting argu¬ 
ments suggest that a PTA with N p pulsars can charac¬ 
terize up to 2N p /7 chirping GW point sources per GW 
frequency bin or 2N p /6 monochromatic sources. This is 
just the number of measurements (an amplitude and a 
phase per pulsar) divided by the number of parameters 
characterizing a single GW source (7 for a chirping binary 
and 6 for a monochromatic binary). We therefore expect 
a PTA to be sensitivity limited when every GW frequency 
bin has more than 2N p /7 sources. At present there are 
more than 20 pulsars in the IPTA with rms timing residu¬ 
als <r rms < l/is, and a few pulsars with cr rms < 100ns [95| . 
With the advent of the Chinese five hundred meter spher¬ 
ical aperture telescope 0 and the Square Kilometer Ar¬ 
ray (SKA) [97}, there will be a major leap in sensitivity. 
A conservative estimate suggests that the SKA could de¬ 
tect more than twenty thousand pulsars, including hun¬ 
dreds of them with <r rms that will match or supersede 
the best pulsars currently known. Such a PTA may no 
longer be a detector capable only of detecting a stochastic 
GW background (i.e., a confusion-limited detector) but 
may become a point source telescope capable of carrying 
out matched-filtering GW searches 94(. In view of this 
bright prospect, we now compute the SNRs of eccentric 
sources in the frequency band of PTAs. 

Since binaries on eccentric orbits radiate in a wide 
spectrum of harmonics n of the mean orbital frequency, 
we can write the SNR as: 


Pt. ( n i /orb) 


hewfo, /orb) 
h 2 ,(n/orb) 


(50) 
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with i = [low, high] and |98j : 


• ( 51 > 

1 + Z 

where h n is given by Eq. ©. Eq. (l5ll) can be inter¬ 
preted as the averaged squared amplitude multiplied by 
the number of cycles completed during the observation 
time T 0 bs- In general, the total SNR of a single resolvable 
source can be written as: 


Umax 

P = 'y ' Plow( n 5 /orb) + y ' Phigh( n ; /orb) , (52) 


where n max is given by n max / orb = /trans- Now, bearing 
in mind that the sensitivity for continuous wave sources 
is given by a piecewise function, let us consider the low 
frequency component. Using Eq. (14711 we find that: 

Umax o 

Plow = C E Q) 9(n,e 0 ) fT > ( 53 ) 

n—l 

C _ 4 n ^ 4/3 N p (N p — 1) rj bg M 1Q / 3 

45 d 2 L {l +z) 2 Ata 2 ms ’ 

(54) 

where we have used f = n f Q ^/(l + z) in the last line. In 
the case where most of the detectable signal is contained 
in modes with n < n max , we can use an analytical form 
for the sum in Eq. (15.3f) . In Appendix lAl we show that: 


°o 2 

G(e 0 ) = E Q) 9(n, e 0 ) 


(1 -elY 


, 85 2 

1 + ¥ 6 ° 


5171 4 , 1751 6 , 297 8 
192 6 ° + l92" e ° + 1024 6 ° 


(55) 


Hence, summing over all the harmonics enables us to 
recast Eq. (El as follows: 


Plow — CG(e$) /o r 6 / 3 ■ (56) 


We can find a similar expression for the high frequency 
contribution, namely: 


„2 ST 9(n,e 0 ) 2/3 

Phigh ° < /n\l '' orb ’ 

„ = „ max («/2) 

- _ 4 y/2^ N p (N p - 1) T obs Ad 10/3 (1 + zf 
45 d 2 L Ata 2 ms 


(57) 


(58) 


If the first harmonic n = 1 is located within the high 
frequency regime (i> /trans), then no detectable signal 
occurs in the low frequency regime, so n max = 1 and 
we can analytically evaluate the sum in Eq. (1571) . In 
Appendix [A] we show that this sum is given by 


y( e o) = E 


n—l 


g(n , e 0 ) 

(n/2) 4 


= 1- -e, 


0 5 


(59) 


and the high frequency contribution can be expressed as: 


Phigh = 6T(e 0 )/ orb 


2/3 


(60) 


We can re-write the low and high frequency contributions 
to the SNR in a convenient way using the transformation 

11 — /orb/ /trans- 


where: 


(1 + z) 2 CG(e 0 )u 16 / 3 , u <C 1, 
(1 + z) 4 CY(eo) u~ 2 ! 3 , u> 1, 


(61) 


_ 4 ^2^/ 3 N p (N p - 1)T^M w / 3 

45 d 2 L Ata 2 ms ■ > 

We note that the requirement that u -C 1 in the first part 
of Eq. (IbTjl is due to the fact that eccentric sources will 
emit in a wide range of harmonics. For more moderate 
eccentricities, this requirement is weakened, such that 
Eq. (1611) applies to all orbital frequencies in the limit of 
very small eccentricity. 

To give a sense of scale, we can reexpress Eq. ED as: 


2 _ ~2 / (1 + z) G(e o) M 16 / 3 , 11 < 1 , 
(I + 2) 4 y(e 0 )u _2/3 , u> 1, 


P =P 


(63) 


p 2 = 4.26 x 10“ 2 N p ( N p - 1) 


T obs \ 5/3 /100 Mpc N 2 


M \ 10/3 

10 8 M q ) 


10 yr/ V d L 
100 ns V ( 0.05 yr 
At 


(64) 


Finally, in the case that individual sources are emitting 
in the transition regime between low and high frequency 
sensitivity (i.e., / or b < /trans, but the eccentricity is large 
enough that significant signal is contained in harmonics 
with n/orb > /trans), the total SNR is given by: 


2 -2 

P =P 


(1 + zY 

N„ 


'‘■max 2 

E (f) g(n,eo)u 16/3 


n—l 


'max / \ 

(i+A E ^“- 2/3 

1 1 (rz/2) 

Umax “I - ! V / / 


(65) 


where formally -/V max —> 00 , but in practice, we find that 
-Nmax = 1500 suffices for all of the eccentricities consid¬ 
ered in this work. In Figure [5] we show the expected 
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SNR p for sources that emit in three different regimes: 
very low frequencies (/ or b ■C /trans) ; transition frequen¬ 
cies (0.01/trans ^ /orb < /trans), and high frequencies 
(/orb ^ /trans)- These results indicate that: 

• Single resolvable binaries that satisfy / or b -C /trans 
undergo a substantial SNR increase. Heuristically, 
we can understand this effect based on the results 
reported in |jj, namely, the SNR gets contributions 
from all harmonics of the orbital frequency n f OI b, 
some of which will be located in the region of max¬ 
imum sensitivity of the PTA. The bottom panel of 
Figure 5 shows that we can analytically compute 
the SNR for binaries with orbital frequencies up 
to /orb <, 0.01/t r ans and eo < 0.8 with an accu¬ 
racy better than 1% using Eqs. (l55l) . (l63l) and (1551) 
for u -C 1. The bottom panel of Figure [5] shows 
that the regime of applicability of these relations 
increases as we consider sources radiating at very 
low frequencies (see the line labelled u = 0.005). 
Using these relations, we find that the increase in 
SNR in the low frequency regime is given by: 


u<0.01 

/^increase 


Pe 0 > 0 
Pe o=0 


\jG(e o) . 


( 66 ) 


Evidently, the contribution from harmonics located 
in the high frequency regime — where the sensitiv¬ 
ity of the PTA is poorer — tends to slow down the 
increase in the SNR and eventually attenuate it. 
This is clearly shown in the top panel of Figure 5. 

• Binaries with orbital frequencies 0.01 /trans 
/orb < /trans need to be described by Eq. (1551) in¬ 
cluding the contribution from harmonics located in 
the low and high sensitivity regime frequency of a 
PTA. In that case we include up to N max = 1500 
to provide a reliable answer. 

• Finally, binaries with / or b > /trans are very well 
described by Eqs. (1551) . (1531) and (1551) for u > 1. 
These relations indicate that the loss in SNR due 
to eccentricity is given by: 


U>1 _ 

^loss 


Pe q>0 
Pe o=0 




(67) 


C. Signal-to-noise ratio calculations for a 
stochastic gravitational wave background 

The nature of a stochastic GW background allows us 
only to predict the statistical properties of the signal it 
generates, not the precise signal. Matched filtering ap¬ 
proaches are not, therefore, applicable and instead we 
rely on cross-correlation of data streams from different 



Eccentricity 


FIG. 5. Expected signaf-to-noise ratio p for sources that may 
be detected in the frequency band of PTAs assuming N p = 10, 
d,L = 100 Mpc, z = 0.022, <T rms = 100 ns, At = 0.05 yr and 
M = 10 9 M 0 (see Eq. (16511 ). The top panel shows the enhance¬ 
ment in p at low frequencies (u -C 1), and the corresponding 
attenuation at higher frequencies. We also compare the per¬ 
formance of the expressions given in Eqs. (ESI) and io) with 
the actual numerical evaluation of Eq. (1651) . 


pulsars. The SNR statistic we shall adopt in this case 
is described in Ref. & This is the linear combination 
of cross-correlations between different pulsars that max¬ 
imizes the SNR, defined as the ratio of the expectation 
value of the statistic in the presence of a signal to the 
rms value in the absence of a signal. The SNR for this 
optimal statistic is 


N v N v 


E2=8 EE T °bs / d f 


i>j 


r 

s 2 M) 


( 68 ) 


For an isotropic background, the overlap reduction func¬ 
tion r,j is entirely determined by the angular separation 
of the pulsars @]. Assuming that the pulsars in the PTA 
are randomly placed on the sky, F^ can be approximated 
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as the rms value over the sky, i.e., 


r« = x = (4v^) 

N p N p 

EE r 


N P (N P - l)x 


*>3 3 

Eq. (1551) thus takes the form 


E 2 = n p (N p -1) r obs | d/ _^(/) 


12 


S*{f )' 


Additionally, 


c 3 H 2 n(f) GW 

Sh{f) = 2n 2 /3 and 

S n (f) = 24n 2 Ata 2 ms f 2 . 


(69) 

(70) 

(71) 


(72) 


1. SNR calculations for binaries with fixed, eccentricity 

The SNR for a population of binaries with fixed orbital 
eccentricity can be derived using Eqs. m, ® and (ED: 

^ _ No 2 M 10 N N p (N p - 1) T obs B 2 (e 0 ) df 
3888 7 t 14 / 3 ( Ata?ms f 4 f 2 */ 3 ' 

l 73 ) 

Using the coordinate transformation v = ///o, with / 0 = 
we obtain: 

iVg .M 10/3 iVp (N p - 1) T Q 2 b 6 s /3 H 2 (e 0 ) 

29808 tt 14 /3 (At af ms ) 2 ‘ 1 ’ 

We thus obtain an expression for the SNR of a stochastic 
GW background of identical constant eccentricities eg: 


/ m \ 10/3 

£ 2 = 23.49 B 2 (e 0 ) A p ( N p - 1) J 

x(^) 26/3 ( N ° 3 ) 2 (— y 

V 10 yr/ \10 -3 Mpc / V ^rms / 


X 


^ 0-05 yr ^ 2 


(75) 


In Figure [D we show the expected SNR from a stochas¬ 
tic GW background generated by sources with fixed to¬ 
tal mass M. These results have been generated using 
the fiducial values quoted in parentheses in Eq. ED, 
and assuming a network of N p = 10 pulsars. Figure [7] 
shows that eccentricity tends to reduce the expected SNR 
from a population of compact binary sources. This ef¬ 
fect is marginal for binaries with low to moderate val¬ 
ues of eccentricity, i.e., for eo € [0, 0.6]. However, the 
expected SNR of a stochastic GW background gener¬ 
ated by a population of highly eccentric binaries satisfies 
E(eo = 0) ;> 10E(e o ~ 0.9). This is a natural conse¬ 
quence of the effect of the attenuation factor B(e o) on 



f A> 


FIG. 6. We show the construction of a function that cap¬ 
tures the harmonic content from fixed eccentricity sources 
with / or b < /o, and that incorporates the contribution from 
evolving eccentricity sources with / orb > /o- The plot shows 
the case for eo = 0.7. Notice that the ‘Physical Evolution’ 
function Z given by Eq. ED reproduces the expected physi¬ 
cal behavior in the appropriate limits. 


the strain of a stochastic GW background (see Figure [2]). 
In the following Section we extend this analysis to con¬ 
sider populations in which the orbital eccentricity of the 
binaries evolves. 



Initial Eccentricity e 0 


FIG. 7. Expected signal-to-noise ratio E for a stochastic grav¬ 
itational wave background generated by sources with total 
mass M = 10 9 Mq and whose eccentricity is either fixed or 
evolving, as indicated in the Figure. We have used the fidu¬ 
cial values quoted in the parentheses of Eq. ED, and assumed 
N p = 10. Note the substantial suppression in signal-to-noise 
ratio due to the effect of eccentricity. As in the case of sin¬ 
gle resolvable sources, eccentricity noticeably suppresses the 
detectability when eo > 0.6. 
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2. SNR calculations for binaries with evolving eccentricity 

A more realistic astrophysical scenario is one in which 
the eccentricity of binaries that generate the stochastic 
GW background is allowed to evolve. In this case, we 
again use Eqs. GD and m, but we now use Eq. m 
along with the function S(f, fo, eo, z ) in Eq. (1M1) to take 
into account the frequency evolution of the eccentric¬ 
ity. Since the function S(f,fo,eo,z) was derived using 
Eq. (l20l) , we identify fo as the orbital frequency at which 
the ensemble of binaries have a fiducial orbital eccentric¬ 
ity e 0 = e(/ 0 rb = /o), where f 0 = T~ h l. 

We can compute the SNR for the evolving eccentricity 
case assuming that the GW background signal evolves 
both above and below / or b = /o, so that e > eo for 
/orb <, fo- In that scenario, the contribution from sources 
for frequencies / or b fo is highly attenuated, as shown 
in Fig. [3] However, this scenario is problematic, partic¬ 
ularly for large values of eo- In reality, we expect some 
dynamical process to be driving binaries to eccentricities 
of eo at /o, so that the behavior of the eccentricity at 
lower frequencies will vary depending on the details of the 
mechanism. In order to make a SNR comparison between 
sources with fixed and evolving eccentricity that does not 
include such severe attenuation for / or b fo, since such 
attenuation is not well astrophysically motivated, we can 
modify the framework described by Eq. (EO). As dis¬ 
cussed by Kocsis and Sesana [39j, the rate of inspiral 
depends on the mechanism driving the evolution, and 
will generically be more rapid than the GW-driven case. 
However, given that the likely dynamical processes pre¬ 
ceding GW domination tend to drive binary eccentrici¬ 
ties to fixed values, one physically reasonable, if simplis¬ 
tic, approach is to assume that sources with / or b < fo 
evolve in frequency at the appropriate rate for gravita¬ 
tional emission, but with constant eccentricity, whereas 
sources with / or b ^ fo evolve following the behavior given 
by the function £(/, fo, eo, z) in Eq. (1M1) . Therefore, the 
attenuation function for this scenario is given by: 


Z(f,fo,e 0 ,z) 


fr/h 

E 


71—1 


1 gi n , e(/o rb ; e 0 )) 
-P(e(/orb; eo)) (n/ 2 ) 2 ^ 


OO 

+ E 

n—fr/So + l 


1 g{n,eo) 
F ( e o) (n/2) 2 ' 3 ' 


(76) 


We show the form of this modified prescription in Fig. © 
assuming a population of sources with eccentricity eo = 
e(/orb = fo) = 0.7. Using this approach, Fig. © shows 
that the expected SNR from sources with evolving ec¬ 
centricity is less attenuated that their fixed eccentricity 
counterparts, which is a natural consequence of the way 
in which we constructed the Z(f, fo, eo, z) function, and 
is the expected physical behavior; since we have found 
that higher eccentricities are more attenuated, the evolv¬ 
ing eccentricity case, which evolves to lower eccentric¬ 
ities due to gravitational-wave emission, should there¬ 


fore be less attenuated than its fixed eccentricity coun¬ 
terpart. Furthermore, evolving eccentricity sources with 
low eccentricities tend to have larger SNR values because 
Z(f, fo,eo,z) ;> 1 for frequencies f r /fo <, 10, and most 
of the SNR is accumulated at lower frequencies due to 
the strong suppression factor f~ 26 / 3 in Eq. (T73l) . Simi¬ 
larly, since highly eccentric systems tend to circularize for 
larger f r /fo values, the net enhancement in SNR of evolv¬ 
ing over fixed eccentricity sources is less pronounced. 

This analysis shows that eccentricity introduces sub¬ 
stantial qualitative and quantitative changes in the prop¬ 
erties of the GWs emitted that, in the context of cur¬ 
rent data analysis algorithms, will make their detection 
more challenging. Developing alternative techniques for 
the detection and characterization of these signals goes 
beyond mere curiosity. Since the orbits of SMBH bina¬ 
ries may only shrink to small enough separations for GW 
domination due to interaction with their environments, 
and these interactions may drive the binaries to large ec¬ 
centricity, it is quite plausible that eccentricity will play 
a fundamental role in the dynamical evolution of SMBH 
binaries within the sensitivity band of PTAs. This arti¬ 
cle is a first step to addressing some of these outstanding 
challenges in the detection of eccentric supermassive bi¬ 
naries. 


V. CONCLUSIONS 

Eccentric binary systems may play a more relevant 
role in the dynamics of compact binary systems than 
previously thought. In light of studies which suggest 
that SMBH binaries may have non negligible eccentricity 
while emitting in the sensitive frequency band of PTAs, 
we have provided a solid foundation to study the prop¬ 
erties of eccentric binary systems. 

In this article we have developed an analytical frame¬ 
work that enables the construction of rapid spectra for 
a stochastic GW background generated by a population 
of eccentric sources which builds upon the work of Phin- 
ney i and Enoki and Nagashima Q. We have also de¬ 
rived several new analytical approximations that expand 
upon the results of Peters and Mathews [l] to fully as¬ 
sess the impact of including eccentricity on the detection 
and characterization of eccentric binary systems in the 
context of single resolvable sources. 

The analytical summations we have derived to bench¬ 
mark the SNR of single binaries that radiate in the high 
frequency regime of PTA sensitivity to continuous wave 
sources do not suffer from the limitations of numeri¬ 
cal summation, particularly for very large eccentricities 
where harmonics at hundreds or thousands of times the 
orbital frequency may significantly contribute to the sig¬ 
nal. Regarding single resolvable binaries that radiate pre¬ 
dominantly in the low frequency PTA sensitivity band, 
our analytical results can be used to benchmark the in¬ 
crease in SNR for sources with eccentricities as high as 
eo ~ 0.8 with an accuracy better than 1%. 




13 


We have provided ready to use expressions to com¬ 
pute the SNR for eccentric single resolvable sources and 
a stochastic GW background generated by a population 
of eccentric binaries. Our results conclusively show that 
eccentricity will have a positive impact on the detec¬ 
tion of single resolvable sources emitting primarily at 
gravitational-wave frequencies / < 2 T7* On the other 
hand, single resolvable sources whose fundamental n = 1 
harmonic is located at a frequency / = / or b > 2T 0 ^, 
or a stochastic, isotropic GW background generated 
by binaries with low to moderate values of eccentricity 
(eo £ [0, 0.6]) may still be recovered with SNRs compa¬ 
rable to their quasi-circular counterparts. The SNRs of 
highly eccentric binaries, however, will be substantially 
suppressed, thus requiring the development of alterna¬ 
tive search techniques to detect and characterize these 
signals. 

In forthcoming work, we will apply the tools developed 
here to devise a new, efficient and accurate framework to 
explore the ability of PTAs to extract the signatures of 
eccentric binary systems and reconstruct the intrinsic pa¬ 
rameters of single resolvable sources and the astrophysi- 
cal distribution of parameters for stochastic signals. 
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Appendix A: Sums of Bessel functions that are 
relevant for the study of eccentric binary systems 

In this Appendix we show how to evaluate the sum over 
all harmonics n for the cases described in the main text 
of the article. The solutions presented in this Appendix 
are based on Bessel’s solu tion of the Kepler equation, 
M = e — e sin E{M, e) [Tonj : 


E(M, e) = M + 2 V Sm(nM) J n (ne ). (Al) 

n=1 

Using the previous relation, we have found the following 
results: 


oo 

E J n( ne ) 

n— 1 
oo 


1 l 

2 + 2(i-e) 1/2 ' 
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n—l ^ (T e ) 
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n—l 
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(A5) 


Using these results and those quoted in the Appendix of 
PM El, we obtain 
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L (e) = E n4 9(n, e) = - - 

n—l (1-e 2 ) 
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(A15) 

(A16) 

(A17) 


Eq. (IA16D was used to derive Eqs. (l^ll- (fTTl) . Eq. (I A 1511 
was used in Eq. ©• The remaining expressions, (IA14I) 
and (I A1 7f) . were used to determine Eqs. (E9l) . (l55l) 

and (1591) . 
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Appendix B: Convergence of infinite sums 

We now estimate how many terms n are needed for 
convergence of sums of the type: 

77 max 

iVKax) = ^ n p g(n, e ). (Bl) 

n—l 

We do this by computing the fractional error in the nu¬ 
merical value of the sum N by including up to rc max har¬ 
monics, and then comparing this value with the exact 
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analytical result and the numerical fit. We consider first 
the well known sum given by Eq. (IA15I) . We have found 
that including up to 100 harmonics is sufficient to repro¬ 
duce the exact analytical result for eccentricities e <; 0.7. 
However, for eccentricities up to e = 0.9 we need to in¬ 
clude up to n = 400 harmonics; n = 800 for eccentricities 
up to e = 0.94 and n = 1200 for eccentricities as high as 
e = 0.96. 

Another important sum is given by Eq. (IA17I) . Fig¬ 
ure [5] shows that this sum is highly convergent. Note 
that n max = 100 harmonics is sufficient to ensure 
that the fractional error <; 0.1% in the entire domain 
e £ [0.0, 0.98]. With n max = 500 harmonics, the 
fractional error <; 0.001% in the entire domain e £ 
[0.0, 0.98]. 




Eccentricity 


FIG. 9. The numerical fit given by Eq. (HO) reproduces 
Eq. (IB 1 1) with p = —2/3 and n max = 1500 with a 
fractional error < 0.01% in the entire domain e £ [0.0, 0.95]. 


FIG. 8. Fractional error in Eq. ED for p — —4. The choice 
n max = 100 is sufficient to ensure that the fractional error < 
0.1% in the entire domain e £ [0.0, 0.95]. 


Appendix C: Attenuation factor B(e) 

We could not find an analytical solution for the fun¬ 
damental sum given in Eq. (|29f) . Instead, we con¬ 
structed a numerical fit that robustly reproduces the sum 
given by Eq. ED with p = —2/3 and n max = 1500 
with a fractional error <J 0.01% in the entire domain 
e £ [0.0, 0.95]. This is shown in Figure [HI 
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